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c3 , Abstract 

We use semiclassical methods to study processes which give rise to change of topology 
and therefore to baryon number violation in the standard model. We consider classically 
allowed processes, i.e. energies above the sphaleron barrier. We develop a computational pro- 
cedure that allows us to solve the Yang Mills equations of motion for spherically symmetric 
configurations and to identify the particle numbers of the in- and out-states. A stochas- 
tic sampling technique is then used to map the region spanned by the topology changing 
solutions in the energy versus incoming particle number plane and, in particular, to deter- 
mine its lower boundary. A lower boundary which approaches small particle number would 
be a strong indication that baryon number violation would occur in high energy collisions, 
Q-i! whereas a lower asymptote at large particle number would be evidence of the contrary. With 
our method and the computational resources we have had at our disposal, we have been able 
to determine the lower boundary up to energies approximately equal to one and a half time 
times the sphaleron energy and observed a 40% decrease in particle number with no sign of 
the particle number leveling off. However encouraging this may be, the decrease in incom- 
ing particle number is only from 50 particles down to approximately 30. Nevertheless, the 
formalism we have established will make it possible to extend the scope of this investigation 
and also to study processes in the classically forbidden region, which we plan to do in the 
future. 
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I. INTRODUCTION 



Since the pioneering work of 't Hooft |l] it has been known that the axial vector 
anomaly implies that baryon number is not conserved in processes which change the topology 
of the gauge fields. Baryon number violating amplitudes are non-perturbative and viable 
methods of calculation are scarce. The two primary methods of obtaining non-perturbative 
information in quantum field theory are either semi-classical techniques or direct lattice 
simulations of the quantum fluctuations. Theories with small coupling constants are not 
suited for the latter, so the electroweak sector of the standard model lies beyond the reach 
of direct lattice calculations. This means that semiclassical methods presently offer the only 
way to study baryon number violating electroweak processes. 

Electroweak baryon number violation is associated with topology change of the gauge 
fields. Classically, gauge field configurations with different topology (i.e. differing by a topo- 
logically non-trivial gauge transformation) are separated by an energy barrier. The (unsta- 
ble) static solution of the classical equations of motion which lies at the top of the energy 
barrier is called the sphaleron 0. At energies lower than the sphaleron energy, topology 
changing transitions, and hence baryon number violation, can only occur via quantum me- 
chanical tunneling. At zero temperature and low energy the tunneling rate can be reliably 
calculated and is exponentially small. A few years ago, however, Ringwald || and Espinosa 
0] noticed that a summation of the semiclassical amplitudes over final states gives rise to 
factors which increase very rapidly with increasing energy. This may lead to a compensation 
of the exponential suppression for energies approaching the energy of the barrier, i.e. the 
sphaleron energy E sp h. Intuitively, one might expect suppression of tunneling to become 
much less severe as the energy approaches the energy of the barrier, in particular, one might 
expect it to disappear altogether for E > -E sp h, i.e. in the region where the topology changing 
processes are classically allowed. Investigations have indeed confirmed that this is precisely 
what happens in high temperature electroweak processes || : as the temperature approaches 
E sp h (which is in fact temperature dependent for a thermal plasma), the barrier-penetration 
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suppression factor becomes progressively less pronounced, and electroweak baryon number 
violation becomes unsuppressed altogether above the critical temperature. The situation 
is, however, much less clear for high energy collisions and it would be premature to con- 
clude that baryon number violation can occur with a non-negligible amplitude. Phase space 
considerations are more subtle and simply because one has enough energy to pass over the 
barrier does not guarantee that one does so. The problem is that in high energy collisions 
the incident state is an exclusive two particle state, which is difficult to incorporate in a 
semiclassical treatment of the transition amplitude. 

A possible remedy to this situation has recently been proposed by Rubakov, Son and 
Tinyakov who suggested that one considers incident coherent states, but constrained so 
that energy and particle number take fixed average values 

E=- 2 (1.1a) 

g 2 

T 

In the limit g — > 0, with e and v held fixed, the path integrals giving the transition amplitudes 
are then dominated by a saddle point configuration which solves the classical equations of 
motion. This permits a semiclassical calculation of the transition rates. Information on high 
energy collision processes with small numbers of incident particles can then be obtained 
from the limit v — + 0. While this limit does not strictly reproduce the exclusive two-particle 
incoming state, under some reasonable assumptions of continuity it can be argued that the 
corresponding transition rates will be equally suppressed or unsuppressed. 

When the energy is below the sphaleron barrier the semiclassical paths that dominate 
the functional integral in Ref. |J must be complex for ( |1.1| ) to be satisfied. Finding such 
solutions is a formidable analytic problem, but one that is well suited to numerical study. 
The numerical evolution naturally divides into two regimes. There is a purely Euclidean 
evolution, corresponding to tunneling under the barrier, and a Minkowski evolution corre- 
sponding to classical motion before and after the tunneling event. The desired semiclassical 
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paths may be obtained by appropriately matching the Euclidean and Minkowski solutions 
onto one another, and the transition amplitude may then be calculated. 

When the energy is greater than the sphaleron barrier, transitions are classically allowed 
and solutions that saturate the functional integral are real. This is the regime examined 
in this paper. When chiral fermions are coupled to gauge and Higgs fields which undergo 
topological transitions, Ref. 0] shows that the anomalous fermion number violation is given 
by the change in Higgs winding number of the classical system. This paper is primarily an 
investigation of whether and to what extent topology change occurs in classical evolution 
with low particle number in the incident state. Since Minkowski evolution is also required 
for the analysis below the sphaleron, the techniques developed in the present investigation 
will be useful there as well. 

The primary impediment for rapid baryon number violation is the phase space mismatch 
between incoming states of low multiplicity and outgoing states of many particles. The 
authors of Ref. || look at simplified models and observe that, classically, it is difficult 
to transfer energy from a small number of hard modes to a large number of soft modes. 
However, the investigations in Ref. |S| find that for pure Yang-Mills theory in 2-dimensions 
the momenta can be dramatically redistributed, although unfortunately the incident particle 
number seems to be rather large in their domain of applicability. Ref. [|l(| studies the Yang- 
Mills-Higgs system in a 2-dimensional plane-wave Ansatz and again finds that momentum 
can be efficiently redistributed. It is the purpose of our investigation to shed further light on 
the situation in 4-dimensions in the presence of a Higgs field and to investigate the relation 
between incoming particle number and topology change. 

Given a typical classical solution, because of the dispersion of the energy, the fields will 
asymptotically approach vacuum values. Consequently, at sufficiently early and late times 
the field equations will reduce to linearized equations describing small oscillations about the 
vacuum and the field evolution will be a superposition of normal mode oscillations. In terms 
of the frequencies u> n and amplitudes a n of these oscillators the energy and particle number 
of ( |1.1|) are given by 
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e = ^2u n \a n \ 2 (1.2a) 

n 

^ = EKI 2 > ( L2b ) 

n 

and we see that for typical classical evolution the energy e and the particle numbers V{ and 
v of the asymptotic incoming and outgoing states are well defined (the energy is of course 
conserved and well defined even in the non-linear regime, although no longer given by ( ]1.2ap ). 
In addition, since the fields approach vacuum values for t — > ±00, the winding numbers of 
incoming and outgoing configurations are also well defined. Because of the sphaleron barrier, 
the energy e of all the classical solutions with a net change of winding number is bounded 
below by the sphaleron energy e sp h- The problem we would like to solve then is whether the 
incoming particle number Vi of these solutions can be arbitrarily small, or more generally, 
we would like to map the region spanned by all possible values of e and z/j for topology 
changing classical evolution. 

One could easily parameterize an initial configuration of the system consisting of in- 
coming waves in the linear regime; however, it would be extremely difficult to adjust the 
parameters to insure that a change of winding number occurs in the course of the subsequent 
evolution. For this reason we will instead parameterize the configuration of the system at 
the moment when a change of topology occurs (this will be our starting configuration), and 
we will then evolve the equations of motion backward in time. Following the time reversed 
evolution until the system reaches the asymptotic linear regime allows us to identify the 
incident particle number z/j. By varying the parameters of the starting configuration with 
a suitable stochastic procedure we will then be able to map the boundary of the region of 
topology changing solutions in the e-v plane. 

Note that the problem of baryon number violation above the barrier may roughly be 
divided into two parts. One must find the set of incoming coherent states which give rise to 
a change in topology of the fields, and one must calculate the overlap between the incident 
two-particle scattering state and such coherent states. Both are very challenging. The 
problem considered in this paper is the more fundamental of the two, in the sense that if 
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topology change cannot occur for coherent states with small average particle number, the 
overlap effect with a two particle beam is a moot point. On the other hand, if a change 
of topology can be induced with arbitrarily low particle number in the incoming state, one 
is at the very least assured that exponential suppression, which is a residual of the barrier 
penetration, will be absent. 

In summary, then, our strategy is the following. We start with a (not necessarily small) 
perturbation about the sphaleron with some energy e. We evolve the configuration until 
it reaches the linear regime, at which time we extract the normal mode amplitudes a n and 
compute the asymptotic particle number v. The time reversed solution will have an incident 
particle number v and will typically undergo topology change, since by construction it will 
pass over the sphaleron barrier. There is of course the possibility that the system will go 
back over the sphaleron barrier and return to the original topological sector, but we check 
against this occurrence by evolving the starting configuration in the opposite direction in 
time and measuring the winding number of the asymptotic state. We can then explore the 
space of topology changing solutions by varying the parameters of the starting configuration 
using suitable stochastic techniques. This permits us to map the allowed e-v plane in an 
attempt to place a reliable lower bound on the incident particle number. If this bound is 
comparable with two particles in the incoming state, it would be an indication that the time 
reversed solution, which passes over the sphaleron barrier, can be excited in a high energy 
collision. Hence, this would be a signal that baryon number violation becomes unsuppressed. 
Likewise, if the bound is large this would indicate that high energy baryon number violation 
is unobservable in a two-particle scattering experiment. 

In what follows we put meat on the bones of the above discussion and present our 
numerical results. The structure of this paper is as follows. In Section II we illustrate the 
general properties of sequences of topology changing field configurations, not necessarily 
solutions to the equations of motion. For simplicity we first consider the 2-dimensional 
Abelian Higgs model. We then examine the 4-dimensional SU(2) Higgs model, but restricted 
to the spherical Ansatz to obtain a computationally tractable system. In Section III we 
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examine the classical evolution in the continuum. Since the field equations are coupled non- 
linear partial differential equations, in Section IV we solve them by numerical techniques. 
In Section V we describe the starting configurations at the moment of topology change, 
i.e. our parameterization of the initial state, and in Section VI we solve the normal mode 
problem necessary for extracting the particle number in the linear regime. In Section VII we 
explain the stochastic sampling technique used to probe the initial configuration space and 
we present our numerical results concerning the region spanned in the e-v plane by topology 
changing solutions. In Section VIII we present concluding remarks and directions for future 
research. The reader who is familiar with the the basic properties of the SU (2) Higgs system 
and of topology changing solutions, and is impatient to learn about our results, may skip 
directly to Section VII. However, in our opinion, much of the value of the research we present 
here is to be found in the formalism we have established to parameterize, evolve and analyze 
classical solutions of the SU (2) Higgs system in the spherical Ansatz. This formalism, which 
is illustrated in Sections II through VI, has not only been crucial for obtaining our current 
results, but we are confident it will be invaluable for further investigation into the problem 
of collision-induced baryon number violation both above and below the sphaleron barrier. 

II. TOPOLOGY CHANGING SEQUENCES OF CONFIGURATIONS 

We start our investigation with the 1+1 dimensional Abelian Higgs system, which 
is defined in terms of a complex scalar field 0(x) and an Abelian gauge potential A^(x) with 
action 



where the indices run over and 1, F^ v = d fl A u — d u A^ and D^cf) = <9 M — iA^cf). We have 
set the coupling constant g = 1 and several inessential constants have been eliminated by a 
suitable choice of units. 

The most important feature of this system is that the vacuum, i.e. the configuration of 
minimum energy, occurs for non- vanishing 0, indeed, in our units for \<p\ = 1. Since this does 
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not specify the phase of 0, there is not a unique vacuum state, but rather multiple vacua. 
Still, because of gauge invariance one must be careful in regard to the physical significance of 
the phase of 0. A local variation in the phase of can always be undone by a suitable gauge 
transformation, and since gauge equivalent configurations must be considered physically 
indistinguishable, local variations of the phase of the scalar field do not lead to different 
vacua. However, variations of the phase of by multiples of 2n (as the coordinate x 1 spans 
the entire spatial axis) cannot be undone by a local gauge transformation, and thus define 
topologically distinct vacuum states. These vacua differ by the global topological properties 
of the field configuration. The condition |0| = 1 restricts the values of the scalar field to the 
unit circle (in the complex plane). In the A = gauge, which we use throughout this paper, 
the values assumed by at x 1 = ±oo stay constant in time. If we demand that takes 
fixed identical values as x 1 — > ±oo (a condition we later relax), then the number of times 
winds around the unit circle as x 1 spans the entire real axis is a topological invariant (the 
winding number) which characterizes different topologically inequivalent vacuum states. 



Imc|) lm<|> lm<|> 




(a) (b) (c) 




(d) (e) (f) 

FIG. 1. Example of two inequivalent vacuum configurations (a, c) and a field configuration at the top 
of the energy barrier separating them (b). Figures a, b and c trace the field in the complex plane as the 
spatial coordinate spans the entire axis. A three dimensional perspective has been added in figures d, e and 
f to illustrate the detailed dependence of on the spatial coordinate. 
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Figures la-c illustrate three possible contours traced in the complex plane by the field 
variable <fi(x l ) as the coordinate x 1 spans the entire space axis. Inequivalent vacuum con- 
figurations with winding numbers and 1 respectively are depicted in Figs, la and lc. In 
the contour of Fig. la the phase of <ft stays fixed at zero as x 1 ranges between — oo and +00, 
whereas it goes once around the unit circle in Fig. lc. Consequently, the corresponding 
vacuum configurations have winding numbers and 1. The detailed variation of the phase 
is immaterial since it can always be changed locally by a gauge transformation. Thus, in 
Fig. la for example, as x varies from —00 to +00 the field does not have to stay fixed, 
but could wander continuously on the unit circle provided the net change in phase is zero. 
However, the configuration of Fig. la cannot be continuously deformed to that of Fig. lc 
without leaving the vacuum manifold. Therefore a continuous path of configurations con- 
necting neighboring vacua must pass over an energy barrier, a configuration which has the 
property that vanishes at a point, rendering its phase there undefined. The smallest such 
energy barrier is called the sphaleron ||, and its Higgs field component is illustrated in 
Fig. lb. Figures ld-f add the additional perspective of spatial dependence for the field 
(/•(a; 1 ). Figures la-c can be viewed as projections onto the complex plane orthogonal to the 
x 1 axis of the curves in Figs. ld-e. 




Rec 




Re<|) 




Re< 



(a) (b) (c) 

FIG. 2. A different gauge equivalent representation of the configurations illustrated in Fig. 1. 

One should note that the periodic boundary conditions on at x 1 = ±00 can be re- 
laxed. Sometimes it is convenient to use the freedom of performing a time independent 
gauge transformation to make 0(oo) and </>(— 00) differ while keeping both fixed in time 
(for solutions, the constancy in time of 0(±oo) follows from the equations of motion in the 
A = gauge). Thus, the configurations of Figs, la-c can be gauge transformed into the 



S 



configurations shown in Figs. 2a-c. In Fig. 2a the phase of <fi changes by — tc as x 1 goes from 
— oo to +00, while in Fig. 2c it rotates by it. As in Fig. 1, the two vacuum configurations 
differ by a phase rotation of 2n, i.e. by a unit change of winding number. In the interme- 
diate configuration (Fig. 2b) the scalar field takes only imaginary values. In this gauge the 
sphaleron configuration takes a very simple form 

(pix 1 ) = i tanh[y/\ (x 1 - c)], = 0, (2.2) 

where c specifies the location of the sphaleron. 

A possible parameterization for the entire evolution illustrated in Fig. 2 can be conve- 
niently written as 

1 + expfir — 2v A (a; 1 — c)] 

A x = ^ , (2 .3b) 

7r cosh[2v A (x 1 — c)j 

with A Q = 0. As the reader can easily verify, for r = —tt/2 and r = 7r/2 the field </> 
reduces to a number of unit modulus precisely spanning the contours of Fig. 2a and Fig. 2c 
respectively (as x l ranges from —00 to +00). The corresponding values of A\ are chosen to 
make the gauge covariant derivative of vanish, thus ensuring vacuum. We should point 
out, however, that ( |2.jj| ) does not represent the solution of any particular set of equations of 
motion (Euclidean or Minkowski). It is merely a compact parameterization of interpolating 
configurations, in terms of two variables c and r, which might be useful in studying sphaleron 
transitions based on the method of collective coordinates. 

Classical solutions of the 2-dimensional Abelian Higgs model can exhibit topology change 
in much the same way as the vacuum-to- vacuum paths described above. If one couples chi- 
ral fermions to the system, the fermionic current has an anomaly which leads to fermion 
number violation in the presence of topology changing classical solutions. Therefore, this 
model would appear to be a very convenient system for a simplified study of baryon num- 
ber violation in high energy processes. However, as we will discuss in a future section, a 
crucial component of the computational investigation is the ability to identify numerically 
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the normal mode amplitudes of the fields in the asymptotic linear regime. No matter how 
non-linear the system may be at any given point in its classical evolution, typically the 
energy will disperse and bring the system to a regime where the fields undergo small os- 
cillations about a vacuum configuration. This dispersion is expected to occur in any field 
theoretical system, unless prevented by conservation laws such as those underlying soliton 
phenomena. Now, while the 2-dimensional Abelian Higgs model does not possess soliton 
solutions, we have observed computationally that the decay of the sphaleron in this system 
nevertheless gives origin to persistent, localized, large oscillations with an extremely small 
damping rate (this observation was also made by Arnold and McLerran in Ref. |T3|). These 
oscillations, illustrated in Fig. 3, make the system quite unwieldy for a computational in- 
vestigation of baryon number violation based on semiclassical techniques. Consequently we 
turn our attention to the more realistic 4-dimensional SU(2) Higgs system. 

__ 3.26525 




FIG. 3. Sphaleron decay in the 2-dimensional Abelian Higgs model: classical evolution of the <fi field. 
The values of the phase of the complex field are coded by shades of gray, and the modulus of the field by 
the height of the surface. The sphaleron decays rather quickly, but leaves behind a quasi-stable oscillating 
remnant. For a full color figure see http://cthulu.bu.edu/~bobs/bviolatc.htm]. 



10 



Throughout this paper we will ignore both the U(l) hypercharge and the back-reaction 
of the fermions on the dynamics of the gauge and Higgs fields. We shall examine the 3+1 
dimensional SU(2) Higgs system, which is defined in terms of a complex doublet $(x) and 
a gauge potential A^x) with action 

S = Jdx* {—Tri^F^ + {D^)W^ - A($ t $ - 1) 2 | , (2.4) 

where the indices run from to 3 and where 

= d^A v - d„A„ - i[A„, A,] (2.5) 
D ll $ = {d ll -iA ll )$ (2.6) 

with A^ = A^a a /2. We use the standard metric rj^ = diag(l, — 1, — 1, — 1), and we have 
eliminated several inessential constants by a suitable choice of units. We have also set the 
coupling constant g — 1, but shall restore it when explicitly needed using the standard 
model value g = 0.652. For our numerical investigation we shall take the Higgs self-coupling 
A = 0.1, which corresponds to mn = 110 GeV. This value of A is small enough that Higgs- 
field dynamics is non-trivial, but large enough to allow many lattice sites to fall within a 
single Higgs Compton wave length. 

Because of the larger dimensionality of space one expects the energy to disperse much 
more readily in this system than in the 1+1 dimensional Abelian Higgs model, an expectation 
borne out by results of Hellmund and Kripfganz ]n| who observed the onset of a linear 



regime following the sphaleron's decay. For a computationally manageable problem, we 
focus on the spherically symmetric configurations of Ratra and Yaffe [pi], which reduce the 



system to an effective 2-dimensional theory. This effective theory, however, still has much in 
common with the full 4-dimensional theory, such as possessing similar topological structure. 
Furthermore, despite its lower dimensionality, we shall see that the effective system still 
linearizes because of explicit kinematic factors of r in the equations of motion (these factors 
are lacking for the 1+1 dimensional Abelian Higgs model). The ease of linearization in this 
effective 2-dimensional theory is physically reasonable since solutions within the spherical 
ansatz can have their energy distributed over expanding spherical shells. 
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Explicitly, the spherical Ansatz is given by expressing the gauge and Higgs fields in terms 
of six real functions a , a± , a , (3 , // and v of r and t: 

A (x,t) = i a (r,t)er ■ x (2.7a) 
A;(x, t) = i [ai(r, t)cr ■ xx 4 + ^M((/ - cr ■ + 1 + ^ ^ e ijk x>a k ] (2.7b) 

2 7" 7^ 

$(x, t) = [/i(r, t) + iis(r, t)<r ■ x]£ , (2.7c) 

where x is the unit three-vector in the radial direction and £ is an arbitrary two-component 
complex unit vector. For the 4-dimensional fields to be regular at the origin, ciq, a, a% — a/r, 
(1 + 0)/r and z/ must vanish like some appropriate power of r as r — > 0. 

Note that configurations in the spherical Ansatz remain in the spherical Ansatz under 
gauge transformations of the form 

A^Ap+itfdJJ /i = 0---3 (2.8) 
$->[f$, (2.9) 

where the gauge function is given by 

U = exp[ifi(r, t)a ■ x/2] . (2.10) 

We require fi(0, t) = to ensure that gauge transformed configurations of regular fields 
remain regular at the origin. This spherical gauge degree of freedom induces a residual U(l) 
gauge invariance in an effective 2-dimensional theory. The action of this effective theory can 
be obtained by inserting (|2.7|) into (|2.4j), from which one finds 



/r°° r 1 
dt J q dr - -r 2 Tf^ + D" X *D^x + r 2 D^*D^ 



2r 2 V J 2 

-Ar 2 (\6\ 2 -l" 



(2.11) 



where the indices now run from to 1 and in contrast to Ref. |14]] are raised and lowered 
with rj^y = diag(l, —1), and where 
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f/w — d^du dvQp (2.12) 

x = a + if3 (2.13) 

<f) = H + iv (2.14) 

DpX = (dp-iaJ X (2.15) 

D,<P=(d,- i -a fl )<f ) . (2.16) 

The action (|2.11|) is indeed invariant under the U(l) gauge transformation 

a^a^ + d^VL (2.17a) 

X^e in X (2.17b) 
0_> e «V 2 0, (2.17c) 

and we see that the spherical Ansatz effectively yields a system very similar to the Abelian 
Higgs model considered above. In this reduced system the variables ao(r, t) and a\(r, t) play 
the role of the 2-dimensional gauge field. The variables x(r, t) and 0(r, t), which parameterize 
the residual components of the 4-dimensional gauge field and the 4-dimensional Higgs field 
respectively, both behave as 2-dimensional Higgs fields. Note that x has a U(l) charge of 
one while <fi has charge one half. Of course, the presence of metric factors (powers of r) in 
the action (|2.11 ) is a reminder that we are really dealing with a 4-dimensional system. 



We shall work in the = (or A = 0) gauge throughout. In the 4-dimensional theory, 
if one compactifies 3-space to S 3 by identifying the points at infinity, it is well known that 
the vacua correspond to the topologically inequivalent ways of mapping S 3 into SU (2) ~ S 3 



15 1 . These maps are characterized by the third homotopy group of SU (2) and a vacuum 



can be labeled by an integer called the homotopy index or winding number. The effective 



2-dimensional theory inherits a corresponding vacuum structure. From ( 2.11 ) it is apparent 
that the vacuum states are characterized by \x\ = \4>\ = 1, with the additional constraint 
that ix*<p 2 = — 1 (as well as D\x = Di<fi = 0). Convenient zero-winding vacua are given 
by Xvac = —i, 0vac = ±1 with ai V ac = 0. There are in fact other vacua with constant fields 
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(and hence zero winding), but from ( |2.7|) they yield singular 4-dimensional fields. Nontrivial 
vacua can be obtained from the trivial vacua via the gauge transformation ( p. 17] ): 

fl/xvac = d^Q (2.18a) 
X™ c =-ie in (2.18b) 
vac =±e^ 2 . (2.18c) 

When 3-space is compactified Q — > 2mr as r — ► oo (for non-zero integers n). Since Q has 
been set to zero at the origin, the winding numbers of such vacua are simply the integers 
n. Note that x vac winds n times around the unit circle while vac only winds by n/2. This 
is because the (f> field has half a unit of U(l) charge while \ has a full unit. Hence, the 
phase change of % is more dramatic in a topological transition, and for this reason we will 
often concentrate our attention upon % rather than 0, even though the Higgs field is more 
fundamental for topology change J?]]. 

As will become apparent shortly, it is often convenient to relax the condition that 3-space 
be compactified. We may then consider vacua ( |2.18j ) for which Q does not become an even 
multiple of ir at large r. In particular, when Q — > (2n + 1)7T, then x vac — > i and </> vac —>■ ±i 
as r — > oo. Then the gauge function U — >• ±icr ■ x and becomes direction dependent, and 
as expected, space cannot be compactified. 

As in the Abelian Higgs model a continuous path in the space of all field configurations 
which interpolates between two inequivalent vacua must necessarily leave the manifold of 
vacuum configurations and pass over an energy barrier. On such a path there will be a con- 
figuration of maximum energy, and of all these maximal energy configurations the sphaleron 
has the lowest energy and represents a saddle point along the energy ridge separating in- 
equivalent vacua 0. In the spherical Ansatzwe can work in a gauge in which the sphaleron 
takes a particularly simple form, with a M = and 

XBph(r)=i[2/(r)-l] (2.19) 
sp h(r) = ih(r) , 
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where / and h vary between and 1 as r changes from to oo and are chosen to minimize 
the energy functional. Note that the field vanishes at the origin and that the \ field 
vanishes at some non-zero value of r. 

This form of the sphaleron, in which the gauge field a M vanishes and the fields \ and 
are pure imaginary, is convenient for numerical calculations. Nevertheless, it is slightly 
peculiar in the following sense. Finite energy configurations, like (|2.19| ), asymptote to pure 
gauge at spatial infinity (note that ^Xs P h0s P h -1 as r oo). Typically a gauge is chosen 
so that the appropriate gauge function is unity at spatial infinity, and then space can be 
compactified to the 3-sphere. But (|2.19| ) gives Xsph — ► i and </> sp h — > i, which as we have 
seen in the discussion following ( |2.18|) corresponds to the direction dependent gauge function 
U — > icr ■ x. So the sphaleron ( |2.19| ) is in a gauge in which 3-space cannot be compactified. 
Note that an arbitrary element of SU(2) can be parameterized by bol + icr ■ b where 1 is 
the two by two unit matrix and b§ + b 2 = 1. Hence SU(2) ~ S 3 , and defining the north 
and south poles by ±1, we see that icr ■ b with b 2 = 1 parameterizes the equatorial sphere. 
Thus the gauge function U maps the sphere at infinity onto the equatorial sphere of SU(2). 
In this gauge, a topology changing transition proceeding over the sphaleron corresponds to 
a transition where the fields wind over the lower hemisphere of SU(2) before the transition 
and over the upper hemisphere after the transition, with a net change in winding number 
still equal to one. The behavior of the x field in a topological transition is then very 
similar to the behavior of the Higgs field in the 2-dimensional model, already illustrated 
in Fig. 2. The behavior of the <j) field is illustrated in Fig. 4. We could of course, and 
sometimes will, work in a gauge consistent with spatial compactification where topological 
transitions interpolate between vacua of definite winding, as in Fig. 1, but the sphaleron 
would look more complicated. The advantage of (|2.19|) from a computational perspective is 
that perturbations about the sphaleron can be more easily parameterized. 
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(a) 



(b) 



(c) 



FIG. 4. Topological transition in the 4-dimensional SU(2) Higgs model: behavior of the <fi field. The x 
field behaves as in Fig. 2. 



III. CLASSICAL EVOLUTION IN THE CONTINUUM 

So far we have only examined topology changing paths that interpolate between 
inequivalent vacua. We are now interested in examining the topological structure of solutions 
to the equations of motion. For vacuum to vacuum sequences it is clear what we mean by 
topology change: this is simply the change in winding number between the initial and final 
vacua. For solutions, however, the situation is not quite so straightforward. Nevertheless, 
topology change can be precisely defined for solutions whose energy density dissipates to 
zero uniformly in the distant past and future, which is the generic case for classical evolution. 
In the asymptotic regime the uniform dissipation of energy renders the system linear and 
the waves can be expressed as small oscillations about vacua of definite winding numbers. 
By the topology change of such a solution, we simply mean the difference in the winding 
number between these two asymptotic vacua. This difference in winding is in fact just given 
by the change in Higgs winding number, and hence is characterized by zeros of the Higgs 
field (although in the spherical Ansatz it is characterized by zeros of both <p and x)- The 
most important physical consequence of this topology change is that when chiral fermions 
are coupled to the system, fermion number violation occurs and is proportional to the change 
in winding of the Higgs field (see Ref. 0). 

We wish to study whether topology change, and hence fermion number violation, can 
occur in the course of classical evolution with small gauge or Higgs particle-number in the 
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incoming state. Since the system we are studying linearizes in the past, the incident particle 
number is defined and our question is well posed. However, the field equations are coupled 
non-linear partial differential equations which we cannot solve in closed form. Our approach, 
then, is to solve the equations numerically with a discretized r-axis and discretized time 
steps, but first it is useful to examine the continuum system. 
The equations of motion obtained from the action ( p.llj ) are 



^(r 2 /^) = i [D,X*X - X*D„x} + ^ r 2 [D v <j>*(j> - <f>*D u 



(3.1a) 



D 2 + M\X\ 2 ~1) + 



1 



X = -; 



(3.1b) 



DV^ + i(| X | 2 + l) + 2Ar 2 (|0| 2 -l 



= i X0* 



(3.1c) 



To solve these equations given an initial configuration, we must specify the appropriate 
boundary conditions. Boundary conditions for the fields at r = can be derived from 
the requirement that the 4-dimensional configurations they parameterize be regular at the 
origin. One finds that the behavior as r — > must be 

a = a ,ir + a ,3^ 3 + • • • (3.2a) 

o-i = ai,o + ai,2 J " 2 + a i,i r + • • • (3.2b) 

a = air + a^r 3 + a$r 5 + . . . (3.2c) 

p = -l+p 2 r 2 + ... (3.2d) 

fi = Ho + [i 2 r 2 + • • • (3.2e) 

v = Vl r + z/ 3 r 3 + . . . , (3.2f) 

where the coefficients of the r-expansion are undetermined functions of time. The r behavior 
of the various fields is determined by the requirement that r = (x 2 + y 2 + z 2 ) 1 / 2 have the 
appropriate power to render 4-dimensional fields analytic in x, y and z. For example, since 
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A is proportional to aocr-x = (a /V)cr-x, ao must be odd in r. In terms of x an d 0, the 
boundary conditions at r = become 

a (0,t) = (3.3a) 

X(0,t) = -i (3.3b) 

Red r <f){0,t) = (3.3c) 

Im0(O,t) = O. (3.3d) 

Since Q(r,t) vanishes at the origin, one can check that these boundary conditions are gauge 
invariant under spherical gauge transformations. 

There is an additional r = boundary condition given by 

Oi, = ai , (3.4) 
which is obtained by requiring that the two terms in ( |2.7b| ) proportional to a ■ x cancel as 



r — > 0. Note that the v = component of ( 3.1a ) is the Gauss' law constraint, and once 



imposed on the initial data it remains satisfied at subsequent times. Substituting (|3.2| ) into 
Gauss' law gives d t (a 10 — a{) = 0. Therefore, if the boundary condition a 10 = aii is satisfied 
by the initial data it remains satisfied. 

We turn now to large-r boundary conditions. Since we are interested in finite energy so- 
lutions, we require that the fields go to pure gauge at large r. Hence, from ( |2.18| ), a M — > d^Q, 
X — > — i exp[iQ] and — ► ± exp[ifi/2] as r — > oo, where Q(r,t) is the spherical gauge func- 
tion defined in ( j2.10| ). We can choose a gauge in which Q at spatial infinity becomes a 
constant, independent of r and t, so that a M — > as r — >• oo. When we compactify 3-space 
and require — > 2n7r at large r for integer n, then x — > — i and 0— >±lasr— >oo. But as 
discussed in the previous section this is inconvenient for parameterizing the sphaleron, and 
instead we will take Q — ► (2n + l)7r for integer n. Then the 4-dimensional gauge function U 
maps spatial infinity onto the equatorial sphere of SU(2), and we cannot compactify space. 
In this case, however, x i an d 0^±iasr^oo. We will choose the plus sign for 0, 
and in summary we take the large-r boundary conditions to be 



a M (r, t) — > (3.5a) 
X(r,t)^i (3.5b) 
<f)(r,t)^i (3.5c) 

as r — ► oo. There will be times in which it is convenient, mostly for purposes of illustration, to 
take the boundary conditions — > 0, \ — >• — i and — > 1 as r — > oo consistent with spatial 
compactification, however, unless otherwise specified, we will use the boundary conditions 

(PD- 

One can now solve the equations of motion for initial configurations and investigate to 
what extent topology changing transitions occur. Since one cannot obtain analytic solutions, 
we will exploit computational methods. These numerical techniques, which are presented 
in the next section, are based on a Hamiltonian formulation, so we close this section with a 
brief exposition of the Hamiltonian approach to the continuous system. 

Central to this approach are the conjugate momenta to the fields, defined by 

1 3C 

E=— — = r 2 (d a 1 -d 1 a ) (3.6a) 

47T OCL\ 

**=hw =D * x (3,6b) 

vr^^WZ^, (3.6c) 



where £ is the Lagrangian density for the action (|2.11 ). Since ao does not appear in (|2.11 ) 



it has no corresponding conjugate momentum and is not considered a dynamical variable. 
Upon inverting ( |3.6| ) for the time derivatives of the dynamical fields, the Hamiltonian of the 
system is found to be H + He where 

E 2 



roo 

H = 4n dr 



o 



2r 



, i + \n x \ 2 + \n <t) \ 2 +\D r x\ 2 + r 2 \D r ^ 2 



+^(lxi 2 -i) 2 + i(lxl 2 + i)^ 2 

+Re(i X *(j> 2 ) + Xr 2 (|0| 2 - I"" 

and 



(3.7a) 
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roo i 

H c = 4:TtJ^ dr a -d r E + - X*^ x ) + ^ W$ ~~ 0* 7r *) 



(3.7b) 



Variation with respect to ao gives Gauss' law 



d r E = z(ir; X ~ X*Kx) + \ K<l> - 0%) • (3-8) 



Note that this is also the v — component of fl3.1a|) . This is not a dynamical equation just 



as cto is not a dynamical variable. In fact, the Hamiltonian formulation makes it clear that 
this equation is a constraint equation and ao is the corresponding Lagrange multiplier. If 
the initial data are chosen to satisfy Gauss' law, it will continue to be satisfied at subsequent 
times. 

In the ao = gauge, the variables 

ai(r), x(r), <P(r) (3.9) 

form a set of canonical coordinates conjugate to the momenta 

E(r) = r 2 <Vi 
*x( r ) = d oX 

^(r) = r 2 <9 o . (3.10) 



The evolution of these variables is generated by the Hamiltonian fl3.7a|) . Gauss' law, Q3.8|) , 
expresses the residual invariance of the system under time independent local gauge trans- 
formations and is imposed as a constraint on the initial configuration. It is subsequently 
conserved by the equations of motion. Given initial data also satisfying the regularity bound- 
ary condition a^o = «i, and using the boundary conditions ( |3.3|) and ( |3.5| ), a regular solution 
is uniquely determined. We now turn to approximating this solution numerically. 
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IV. CLASSICAL EVOLUTION ON THE LATTICE 



To solve the equations of motion numerically the system must be discretized. For 
this purpose we subdivide the r-axis into N equal subintervals of length Ar with finite 
length L = N Ar. Thus, the lattice sites have spatial coordinates ri = iAr with % = • • • N 
(for our numerical simulations we shall take N = 2239 and Ar = 0.04, giving a lattice of 
size L = 89.56). It is convenient to use the formalism of lattice gauge theories in assigning 
the space components of the gauge fields to the oriented links between neighboring sites 
and in the definition of gauge-covariant finite difference operators. For simplicity, we will 
identify the lattice links via the midpoints between lattice sites, which have coordinates 
r i+1/2 = (i+ 1/2) Ar with i = 0---JV-l. 

The variables for the discretized system will now be defined as follows. The zero com- 
ponent gauge degrees of freedom are defined over the lattice sites, and are given by 

ao,i(t) for i = l...JV-l (4.1) 

with ao,o = clo,n = 0. The spatial components of the gauge field are defined over the links 
of the lattice. We will use the notation aij, or simply a^, to represent the gauge variable 
defined over the link between r\j and r i+1 . This gives the variables 

Oi(t) = a hi (t) for i = 0...JV-l. (4.2) 

As we show momentarily, boundary conditions for the spatial variables a« are not required 
to determine the evolution of the system. However, just as in the continuum, we will impose 



an initial data boundary condition on a corresponding to (|3.4 ) to ensure the regularity of 
the four dimensional fields at the origin (this condition will be discussed shortly). 
The other field variables become 

Xi(t) for i = l...JV-l. (4.3) 

with xo = — i, Xn = i and 
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i(t) for i = l...JV-l (4.4) 



with (fix = i- We are using boundary conditions at r = L motivated by ( |3.5| ). These bound- 
ary conditions do not admit spatial compactification and are chosen so that perturbations 
about the sphaleron may be parameterized more conveniently. Occasionally we will take the 
boundary conditions xn — ~i and 0^ = 1 consistent with spatial compactification; however, 
unless otherwise specified we will use the aforementioned large-r boundary conditions. 

The value of <fi at r = has so far not been specified. We will return to this in a 
moment, but first we consider the discretized covariant derivative. The time-like covariant 
derivatives need no modification, but the continuum covariant spatial derivatives are replaced 
by covariant finite differences, e.g. 

expl — ia; Arl Yj+i — Yi , 
D rX -^— A ^ — — i = 0---iV-l, (4.5) 

and like the gauge fields they are to be thought of as being defined on the links between 
lattice sites. The rest of the discretization is straightforward, and one obtains a discretized 
action Sd expressed in terms of a finite set of variables which still possess an exact local 
gauge invariance: 

a>o,i -> o-o,i + <9A i — 0...N (4.6a) 

a i ^a i + {Q t+1 -Q l )/Ar i = . . . N — 1 (4.6b) 
Xi^e^Xi i = 0...N (4.6c) 

<f>i -> e^ /2 (j) t i = . . . N . (4.6d) 

The discretized gauge function f^(t) with i — . . . N is defined over the lattice sites, and 
satisfies Qo(t) = to maintain the regularity of the corresponding 4-dimensional gauge 
transformed fields. 

Before we continue, however, we must derive the boundary condition for </> at % — 0. 
This is obtained from ( |3.3c| ) and ( |3.3d| ), in which the continuum field at r = is real with 



vanishing spatial derivative. Since a statement about the "derivative" is not gauge covariant, 
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we prefer to state that the real part of the covariant derivative <9 r — iacj), together with the 
imaginary part of 0, must vanish at r = 0. This is equivalent to ( |3.3cj ) and ( |3.3d| ) since 



is real at r = 0. But it has the advantage that it translates into the following boundary 
conditions for the discretized case: 

Re[exp( ~' Q 2 ° Ar )0 1 - 0o ] = (4.7a) 

Im0 o = , (4.7b) 

where a® is the value of a\ t i at % = and should not be confused with the time-like vector 
field. Thus, we write the boundary condition as 

„ r ,—icLn Ar , . 
0o = Re[exp( £ )0 X ] , (4.8) 



which allows us to eliminate 0o from the list of dynamical variables. 
The discretized Lagrangian becomes 

L = 47T £ (dM - a °^-^ V _ |expHa,Ar) X , +1 -X,| 2 l Ar 

~^ I 2 \ Ar / Ar 2 J 

N-l , ■ I ■ A„/n\A A 12 



+ 4tt J] ||(^-ia 0)i ) Xi |2 + r 3|(^- 

i=l 



ia ^ 2 2 | exp(-zaj Ar/2)0 i+ i - 0^ 



2~^<| -r m/2 ^ 

~ (IX.I 2 + 1)|0,| 2 - Re( % *0 2 ) - ^(N 2 - I) 2 - Ar 2 (|0,| 2 - l) 2 }Ar 
_ 47rr 2 [Im(exp(-za o Ar/2)0 1 )] 2 ^ 

This Lagrangian was obtained by discretizing the system as previously explained and by 
replacing 0o by the right hand side of ( |4.8|) . One might think this induces an additional 
contribution to the kinetic term of 0i from the time derivative of ( |4.8|) . However, the term 
proportional to 0o vanishes since it is multiplied by Tq = 0, and hence ( |4.9| ) is the complete 
Lagrangian. 

We define conjugate momenta (the factor l/47rAr is introduced so as to have Poisson 
brackets with a continuum like normalization { r K*,4>j} = Sij/Ar etc.) 
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TP dL _ „2 /p. a 0,»+l «0, 

AnAr d(o ai) ' Ar 



1 dL 
p i = = o 

4nAr d(d a / 



1 9L A 

4vrAr <9(<9 x*) 



i = . 


.N — 1 


(4.10a) 


i = 0. 


. iV 


(4.10b) 


i = 0. 


.AT 


(4.10c) 


i = 0. 


.A^ . 


(4.10d) 



Equation (|4.10b|) is a primary constraint equation, in the sense of Dirac. From (|4.9|) and 
Q4.10Q we obtain the Hamiltonian H + He, with 

Ef |exp(- ! a,Ar) x , +1 -x, l2 > 



i=0 ^'1+1/2 ^' J i=l 



|2 



7TV 



' 1 



2 |exp(-ia 4 Ar/2)0 i+1 - 0i| 2 1 2 2 2 

+^+1/2 ^2 +2(1^1 + R W A 



IT l J 

2 



- l) 2 + Ar?(|tf ~ D 3 } Ar + ^r\ n MfgH^MWl! (4 . lla) 



and 

iv-i ( E — E % 1 

tf c = 4vr J] ao,i{ - * Ar — + i(p?X< " + 5 - #*0} A ^ • ( 4 - llb ) 

Upon commuting (or more precisely, taking the Poisson bracket) the constraint ( |4.10bjj with 
H + He one obtains as a further constraint Gauss' law 

El ~ Et ~ l = i(p* Xt - X*iPi) + l M<f>i - fin) = ji j i = 1 • ■ • N - 1 . (4.12) 
Ar 2 

We impose the second class constraint ao,i = for i = 1 . . . N — 1. The equations of evolution 
that follow from H are then 

TP 

i = 0...N-l (4.13a) 

(4.13b) 
(4.13c) 

and 



ddi 






~dt 


T i+l/2 




dXi 






dt 


= Pi 


i 


dtfii 


_ VTi 




dt 


' i 
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dEi _ exp(zo< Ar)x-i - xl exp(-ia; Ar)x l+ i 



dt 



dpi 
dt 



dt 



+ i 



Ar 

r i+i/2 0*+i exp(i ai Ar/2)0i - <j>* exp(-ia t Ar/2)0 m 



exp(-ia,i Ar)xi+\ 



Xi 



Ar 2 

Xi\<pi\ 2 + i<Pi 



+ 



Ar 

exp(iai_i Ar)xi-i - Xi 



0...iV 



Ar 2 



r 2 

'i+1/2 



exp(-mj Ar/2)0 i+ i 



iV 



Ar 2 



+ r, 



i— 1/2 " 



exp(zai_i Ar/2)0j_! - ^ 



Ar 2 



<MN 2 + i) 



+ i X i(j>*-2\r 2 i (j) i {\(j) i 



r 



N - 1 



(4.14a) 



(4.14b) 



(4.14c) 



where O is given by (|4.8|) , 4>n = h Xo — ~i an d Xn = « (or — — i an d 0w = 1> if as we 
will occasionally do, boundary conditions consistent with spatial compactification are used). 
The momenta of x an d vanish at i = and i = N. 

In summary, we have the following table of independent dynamical variables and their 
respective conjugate momenta: 



variable 


momentum 


index range 


number 


at 


Ei 


i = 0...JV-l 


N 


Xi 


Pi 


i= 1...JV- 1 


2{N- 1) 


<j>% 




i= 1...JV-1 


2(JV- 1) 



Since we have set a 0) j to zero, the number of dynamical variables and momenta (excluding 
boundary fields at r = and r = L) are 2(5iV — 4). Note that ( |4. 13| ) and (|4.14j) give 
2(5A^ — 4) equations, so the system is uniquely determined given the initial values of the 
fields x an d and their momenta (note that boundary conditions for the spatial gauge field 



a,i are not required). The initial data must be chosen to be consistent with Gauss' law ( 4.12| ). 
We will also impose the boundary condition a = Re(xi — Xo)/ 1 Ar, which approximates the 
continuum relation ( |3.4j ). (This relation, which would be conserved in the continuum limit, 
will remain satisfied to O(Ar) in the evolution of the discretized system). 

The restriction to uniform spacing of the subintervals on the r-axis is not fundamental 
and we have also implemented a discretization in which Ar increases as one moves out on 
the r-axis. In this manner one can effectively make the system larger and delay the effects 
of the impact of the waves with the boundary without worsening the spatial resolution near 
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r = 0, where most of the non-linear dynamics takes place. We have found, however, that 
the advantages one gains hardly warrant the additional complications introduced by the 
non-uniform spacing. 

For the numerical integration of the time evolution we have used the leap-frog algorithm. 
Since this algorithm constitutes one of the fundamental techniques for the integration of 
ordinary differential equations of the Hamiltonian type and as such is textbook material, we 
will not discuss it in depth. Essentially, given conjugate canonical variables and pi which 
obey equations 

dqi 

d f t = m , (4.15) 

one evolves the values of q and p from some initial t to t + At as follows. In a first step Pi 
is evolved to the mid-point of the time interval by 

Pi ->p'i =Pi + fi(q)^Y 

= * (4.16) 

(although qi is left unchanged, it is convenient to consider the step formally as a transfor- 
mation of the entire set of canonical variables). In a second step one evolves the coordinates 
from their initial value qi = q[ to their value at the end of the interval 

p'i^Pi = Pi 

4^<g = <fi + 9iW)*t. (4.17) 
Finally, the momenta are evolved from their value at the midpoint to the final value 

pI^pT-p'I + W)^ 

41 - gf = • (4-18) 

One can easily verify that these equations reproduce the correct continuum evolution 
from t to t + At up to errors of order (At) 3 . Moreover, the algorithm has the very nice 
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property that all three steps above constitute a canonical transformation and that it is 
reversible (in the sense that starting from q'", —pf, up to round-off errors one would end 
up exactly with qi, —pi). Because the physical solutions of interest are the time reversed 
processes of the ones we numerically evolve, it is important that we use an algorithm that 
is reversible. Another very nice feature of the algorithm is that, although the evolution of 
the variables is affected by errors of order (At) 3 , the energy of a harmonic oscillator, and 
therefore of any system which can be decomposed into a linear superposition of harmonic 
oscillators, is conserved exactly (always up to round-off errors, but if one works as we do 
in double precision, these are very small). Since extracting the asymptotic normal mode 
amplitudes is the heart of our numerical approach, it is also important to have an algorithm 
that is well behaved in the linear regime. One final comment is in order. In a sequence 
of several iterations of the algorithm, after the momenta have been evolved by the initial 
At/2, the first and third steps, Q4.16|) and ( |4.18|) respectively, can be combined into a single 



step, whereby the momenta are evolved from the midpoint of one interval to the midpoint of 
the next one "hopping over" the coordinates, which are evolved from endpoint to endpoint. 
This motivates the name assigned to the algorithm. 

V. THE INITIAL CONFIGURATION: PERTURBATION ABOUT THE 

SPHALERON 

With a good grasp on numerical solutions of the equations of motion, we can 
turn now to the second crucial component of the computation, namely the parameterization 
of the initial configuration. One could easily construct an initial state consisting of an 
incoming wave in the linear regime; however, it would be very difficult to ensure that such 
a configuration underwent a topology change during its subsequent evolution. Instead, it is 
much more convenient to parameterize the initial state at or near the instant of topology 
change. The system is then allowed to evolve until the linear regime is reached, at which 
point the particle number can be extracted in the manner explained in the next section. 
The physical process of interest is then the time reversed solution, which starts in the linear 
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regime with a known particle number and undergoes a change of topology at subsequent 
times. (In fact, it must be explicitly checked that the winding number of the outgoing 
configuration is different from the incoming one, ensuring that the topology has changed, 
since the system could pass back over the sphaleron barrier and into the original topological 
sector. We have found however that topology change does typically occur.) 

Topology changing transitions within the spherical Ansatz are characterized by the van- 
ishing of at r = and the vanishing of x & t nonzero r. The zero of \ is reminiscent of the 
zero which characterizes the sphaleron of the Abelian Higgs model. However, as shown in 
Ref. 0, it is the zero of the Higgs field (i.e. the zero of <fr) which carries a deeper significance 
and should be associated with the actual occurrence of the topological transition. For a 
sequence of configurations that pass directly through the sphaleron these two zeros occur 
at the same time. Nonetheless, this is not the most general case and the zeros of <fi and \ 
need not occur simultaneously (although for a topological transition, both fields will vanish 
sometime during their evolution) |Tj|. We are free then to parameterize initial topology 
changing configurations imposing that either cf) vanish at the origin or that \ has a zero 
at some non-zero r. It is convenient to choose the latter, in which we parameterize the 
initial configuration in terms of coefficients c n of some suitable expansion of the fields and 
their conjugate momenta, constrained only by the boundary conditions and the requirement 
that the field x has a zero at some non-zero r. Furthermore, we can use the residual time 
independent gauge invariance to make x pure imaginary at the initial time. The field 
is only restricted to obey the boundary conditions and does not necessarily vanish at the 
origin (although it will vanish at the origin at some instant in its evolution if the topology 
is to change). 

To be more specific, we parameterize each field as a (not necessarily small) perturbation 
about the sphaleron given by a linear combination of spherical Bessel functions with the 
appropriate small-r behavior of ( |3.2| ). We only need the first three functions, 




(5.1a) 



x 




(5.1b) 



(5.1c) 



28 



since jo(x) ~ 1, ji(x) ~ x and j2{x) ~ x 2 at small x. Motivated by the boundary conditions 
( |3.5|) , we require the perturbation to vanish at r — L. We thus parameterize perturbations 
about the sphaleron in terms of j nm (r) = j n (a nm r / L) with n = 0, 1 or 2, where a nm are the 
zeros of j n (x), i.e. j„(a nm ) = with m — 1, 2, • • •. The functions j nm (r) form a complete set 
for every n, and the small-r behavior determines the appropriate value of n for each field. 
The reader should note that the expansion of the starting configuration in terms of Bessel 
functions is largely a matter of convenience. This expansion is not related to the expansion 
of the fields in the linear regime (to be discussed in the next section), and any complete set 
of functions with the correct behavior as r — > can be used to parameterize a perturbation 
of the sphaleron localized in the neighborhood of the origin. 

Recall that we must impose the boundary condition a 10 = a\ on the initial data (using 
continuum notation). We are working in the a$ = gauge, but we still have the freedom to 
impose a time independent gauge transformation on the starting configuration to set a = 0. 
Therefore, ( |3.2b[ ) gives ai(r) ~ r 2 at small r, and hence ai(r) is expanded only in terms of 



j2{x). We are thus led to parameterize the initial configuration by 

A^sph 

X(r) = Xsph(r) c i™ hm(r) 

-^sph -^sph 

0(r) = sph (r) + c 2m jom(r) + i c 3m ji m (r) 

m=l m=l 

-^sph -^sph 

Q>m J2m(^) 

m=l m=l 
N i N u 

Y C 6"> j0m( r ) +i Y 

^7m Jlm(^) 

m=l m=l 

fl i( r ) = Yl c » m hm{r) , 



(5.2a) 
(5.2b) 
(5.2c) 
(5.2d) 
(5.2e) 



m=l 



where Xsph — *(2/ — 1) and sp h = ih as in ( |2.19D , and where we have cut off the sums 



at some iV S ph < N. The most general initial configuration is obtained with A^ sp h = N, 
but to avoid exciting short wavelength modes which only correspond to lattice artifacts, we 
take iV S ph < iV/5 to N/10. This implies no limitations on the physical properties of the 

29 



system other that those coming from an ultraviolet cutoff (finite Ar) anyway, and as one 
expects this is born out by numerical results in which typical solutions excite only modes 
with wavelength substantially larger than the lattice spacing. As the dimension of the initial 
configuration space is 8-/V S ph, and since the lattice we work with is rather large, to improve 
the efficiency of our stochastic search we have taken iV S ph ~ iV/50 (iV S ph = 50 for N = 2239). 

To obtain the correct small-r behavior of 7r^, we have inserted an explicit factor of r 2 in 
fl5.2d|) because = r 2 d (p. The profile functions / and h satisfy the boundary conditions 
/(0) = h(0) = and f(L) = h(L) = 1, and will be specified momentarily. For now it is 
sufficient to note that since x(0) = — % and x{L) = i, and since x{ r ) is pure imaginary, it will 
necessarily have a zero for some r > 0. Hence, (|5.2|) specifies a configuration at the moment 
in which \ vanishes. We should also point out that because of its large-r behavior, (|5.2|) is 
expressed in a gauge in that is inconsistent with spatial compactification. 

We have so far used continuum notation, but (|5.2j ) is to be understood as determining 
the configuration at the lattice sites r = r 4 for (|5.2a| )-( |5".2d| ) and at r = r i+1 / 2 for ( |5.2ej ), 

i.e. Xi = x( r i), Pi = 7Tx( r i)) 4>t = 4>( r i), n = ^4>( r i) and a i = ai(r i+1/2 ). We have not yet 
specified the electric field, but since the initial configuration must satisfy Gauss' law we can 
determine by integrating ( |4.12|) outward from i = §toi = N — 1. The value of E must 



be given for this procedure however. In the continuum E(r = 0) = 0, so one is tempted to 
set Eq — 0. But since Eq lives on the first link at r = 7*1/2 = Ar/2, it is better to set 

Ar jp+ji _ Ar [i(p{xi ~ X1P1) + ^K0i ~ <Pl^i)/ 2 } ( r 3) 

2 2 4 > ' 

and then subsequent values of E^ for i — 1- ■ ■ N — 1 can be obtained by integrating ( |4.12|) . 



The sphaleron Xs P h, sp h of (|2.19|) is parameterized by profile functions f(r) and h(r) 
and is a saddle point of the potential energy functional with one unstable direction. This 
direction involves an excitation of the 2-dimensional gauge potential a\. Hence the sphaleron 
is an absolute minimum of the potential obtained from ( |4. 1 la| ) by dropping the a% terms 
(and all the momenta). Using the method of conjugate gradients, with an initial guess for 
/ and h that satisfies the appropriate boundary conditions, we can obtain an extremely 
accurate approximation to the sphaleron by minimizing 
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N ~ 1 f Iv v 12 

^sph/47r = 2^ { — \-r 4 



i=0 



Ar 2 



i+i/2 Ar2 1 + 3(1*^+ w+ 

. JV-l x 

Re(ixJ#) + A r 2 (I^l 2 - l) 2 Ar + £ -^(| Xi | 2 - l) 2 Ar 

J i=l Zr i 



N-l 



g { 4(/ ' A ; /,)2 + -t^^^ + 2 « - w + 



AT-l 



+Ar 2 (/r 2 -l) 2 Ar+E^(l-/,) 2 A 



i=i r * 



(5.4a) 



(5.4b) 



where we have used the boundary condition <fio = Im0i to extend the sum on in (|4.11a|) 
to include % = 0. In our units and with <7 = 1, the energy of the sphaleron is then given by 
e sph /47r = 2.5426 for A = 0.1. 




FIG. 5. Decay of a small perturbation about the sphaleron: evolution of the x held. The values of the 
phase of the complex field are coded by different shades of gray, and the modulus of the field by the height 
of the surface. As explained in the next section, the asymptotic linear system has a particle number of order 
53. The lattice parameters are N = 2239 and Ar = 0.04 with a Higgs coupling of A = 0.1. The initial 
configuration is given by (5.2) with c^i = 0.00247 being the only nonzero parameter. For a full color figure 



see http: / /cthulu. bu.edu/~bobs /bviolatc.htm] 



We are now in a position to numerically evolve perturbations about the sphaleron. Fig- 
ure 5 illustrates the behavior of the x field for an initial configuration given by Q5.2| ) with 
C4 m= i = 0.00247 and all other c-parameters zero. This is in fact the configuration from 
which we have chosen to seed the stochastic sampling procedure which we will describe in 
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Sec. VII. We have found it very convenient and informative to use color to code the phase 
of the complex fields. Unfortunately the illustrations in these pages cannot be reproduced 
in color and we have tried to render the variation of the phase with a gray scale. At some 
point a gauge transformation has been performed in Fig. 5 bringing the asymptotic lin- 
ear state into the sector of zero winding number (consistent with spatial compactification). 
The gauge transformation is made manifest by the sudden change of shading of the sur- 
face. We have performed this gauge transformation because eventually we want to study 
the topology change of the time reversed solution (cf. Fig. 6 below), and this is best done 
in a gauge in which the asymptotic linear state has zero winding number. Moreover, the 
gauge transformation also serves to give a graphic illustration of the gauge invariance of our 
procedure, which is made manifest by the fact that although the shading (or color) of the 
surface changes, there is no discontinuity in the surface itself. 

From Fig. 5 it is clear that the energy, which is concentrated in the neighborhood of 
r = 0, disperses and gives rise to a pattern of outgoing waves. The waves soon become 
linear and possess a definite particle number, in this case of order 53 physical particles 
(using units appropriate to the standard model, which we will refer to as physical units). 

The physical process of interest is then the time reversed solution which starts in the 
linear regime with known particle number, proceeds through the non-linear sphaleron per- 
turbation (|5.2|) at intermediate times and finally linearizes once again at late times. Because 
of time invariance of the equations of motion, this process can be obtained by first evolving 
the perturbation ( |5.2| ) until the linear regime is reached, and then reversing the momenta 
and evolving that configuration forward in time. The resulting solution retraces the evo- 
lution of the sphaleron decay, and then proceeds over the barrier into another topological 
sector. Since our numerical strategy for obtaining asymptotically linear topology changing 



solutions relies upon first evolving the sphaleron perturbation, we shall refer to ( |5.2j ) as the 
"initial" state, while the asymptotic linear states of the physical process will be called the 
"in" and "out" states. 
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FIG. 6. Topology changing transition: behavior of the x field obtained from Fig. 5 by the time reversal 
procedure described in the text. The various shades of gray code the phase of the complex field. The field 
starts as an excitation about the trivial vacuum, passes over the sphaleron and then emerges as an excitation 
about the vacuum of unit winding. Note the persistent strip of 2tt phase change near r = after the wave 



bounces off the origin. For a full color figure see http: / /cthulu. bu.edu / ~bobs /bviolatc.htm] 



Figure 6 represents a physical process obtained from Fig. 5 in the above manner, and it 
illustrates the evolution of the x field for a topology changing solution. The "initial" state 
in Fig. 6, determined from ( |5.2D by the coefficients c n , corresponds to the time-slice half 
way through the depicted evolution. We have reverted to a gauge in which the boundary 
conditions are xn = —i and <pN = 1, consistent with spatial compactification, and in 
which the in-state has no winding and the out-state has unit winding number. This process 
represents an imploding spherical energy shell that converges on the origin, where a change 
of topology takes place. The topology change is indicated by the strip of rapidly varying 
tonality which persists in the neighborhood of the origin and codes the variation of the 2-7T 
phase change of x- With color, this strip would appear as a vivid rainbow, left over as a 
marker of the change of topology of the evolving fields. 

It is important to keep in mind that an arbitrary configuration ( |5.2|) does not necessarily 
produce a topology changing solution, in the sense that at late times the out-state might 



evolve back into the original topological sector. With our parameterization (|5.2| ), however, 
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we have found that the system does in fact typically change topology. Nonetheless, using 
the time reversed procedure above we can always verify whether the in- and out-states have 
the same topology, and if so the initial configuration that produced them can be rejected (or 
equivalently, and more efficiently, we can evolve the initial configuration ( |5.2| ) both forward 
and backward in time and compare the asymptotic states obtained in this way). 

We now have a procedure for constructing solutions which, in the course of their evolu- 
tion, undergo changes of topology. By varying the values of the parameters c n we will be 
able to study the properties of such field evolution and, in particular, explore the domain of 
permissible values for e and v. Before we can implement this procedure, however, we must 
devise a way to calculate the particle number in the asymptotic linear regime. In the next 
section we describe how this can be done. 

VI. NORMAL MODES 

Given an initial configuration parameterized by the coefficients c n , we evolve the 
system until the linear regime is reached, where the fields undergo small oscillations about 
a vacuum configuration. The normal mode amplitudes a n may then be extracted and the 
particle number computed using (|1.2b|) . We turn now to the problem of identifying the 



normal modes. 

Since we have put the field theoretic system of interest on a spatial lattice, to be entirely 
consistent we should also solve the normal modes problem on the lattice. The discrete 
problem, however, cannot be solved analytically and one must resort to numerical methods. 
On the other hand, the normal modes of the continuum system, even restricted to a box of 
finite size L = NAr, can be found analytically. We have solved the problem both numerically 
on the lattice and analytically in the continuum limit. The lattice we consider (N = 2239 
with Ar = 0.04) is big enough that there is excellent numerical agreement between the 
normal modes found by the two methods (the difference between the normalized modes 
never exceeds 10~ 6 ), so we will present here only the continuum solution. 
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Following Ref . [|L6] , we work in terms of gauge invariant variables. We write the fields \ 



and <p in polar form, 

X = -*pe w (6.1a) 
<P = ae iri . (6.1b) 

The variables p and a are gauge invariant. We can also define the gauge invariant angle 

£ = 6-21). (6.2) 

Finally, in 1+1 dimensions we can write 

r 2 /^ = (6.3) 

where eoi = +1 and p, v run over and 1. The variable ip is gauge invariant. Rather 
than working with the six gauge-variant degrees of freedom x, <p and a M we use the 



four gauge invariant variables p, a, ip an d £• These variables satisfy the equations [16 



O^p - ^ 4 {p 4\ r2a 7y + -i(P 2 ~ ^)P+- 2 P° 2 - 2 ^ cos£ = (6.4a) 
d y d , a _ \r 2 a{^d^-^^f + l_ {p2 + + 2Ar2 ^ _ x y _ pacos ^ = Q (g 4b) 



t P + V ° J r 



a , ] p 2 {\r 2 ° 2 d^ + e, u d^) \ 1 2 



where the indices run over and 1 and are raised and lowered with the metric@ 
rjfja, = diag(l, —1), so that c^c^ = df — <9 2 . The energy takes the form 



1 The sign convention of the metric in this paper is opposite to that of Ref. 
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e = 4^°° dr{ (<V) 2 + (<9 r p) 2 + r 2 {d t af + r 2 {d r af + 



12 2 2 

^r^cHp" 1 

p2 _|_ I r 2 a 2 



2^ 2 



(ftO 2 + (drO 2 
1 



+ 



p + 4r 2 cx 2 



+ 



+ ^(P 2 ~ 1) + 2<> + ^ " P ff cos£ + Ar 2 - 1 



(6.5) 



and we see that the vacuum is given by p vac = 1, cr vac = 1, = and £ vac = 0. 

We wish to consider small fluctuations about the vacuum. It is convenient to define 
shifted fields y and h by 



p(r,t) = l + y(r,t) 

ffM) = 1 + MM) 



Then to linear order in h, y, ip and £, 



becomes 



d u &* + AX )h = 



l + ±r 2 



-6 = 



(6.6a) 
(6.6b) 



(6.7a) 
(6.7b) 
(6.7c) 
(6.7d) 



Equation ( |6.7a| ) corresponds to a pure Higgs field excitation characterized by mass 
ran = 2V% while (|6.7b|) -( p.7cl| ) are the three gauge modes of mass m w = l/\/2.i To 
implement the boundary conditions (|3.3j ), we take the gauge invariant fields h, y, ip and £ 
to vanish at r = 0. At r = L we take h, y and £ to vanish (consistent with x an d <fi taking 
their vacuum values there). The r = L boundary condition on ip is that d r ip is zero (ip can 
not vanish at large r since it is proportional to the time derivative of the gauge field). We 



2 Upon restoring the factors of g and the Higgs field vacuum expectation value v, these masses 
take the standard form m# = V2Xv and mw = 
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wish to solve (|6.7|) subject to these boundary conditions, and then extract the corresponding 
amplitudes. 

Let us examine the four types of modes in turn. They can all be expressed in terms 
of the spherical Bessel functions ( |5.1| ). Equation ( |6.7a| ) produces an eigenmode whose non- 
vanishing components are of the form h n (r,t) = h n {r) cosc<Ji n t, with 



h n {r) = Ai„r j (Xi n r) N ln , 



(6.8) 



where uj\ n = (AX + X\ n ) 1 ^ 2 and Ai„ = nir/L for n — 1, 2, • • -. The parameters X\ n have been 
chosen so that h n (L,t) = 0, and the normalization constants N\ n are taken to be 



N- 



2 
I 



1/2 



(6.9) 



so that the h n (r) are orthonormal over the interval [0, L\. To extract these modes from a 
given solution we expand the Higgs excitation as 



h(r, t) = ^2 A n h n (r) cosu ln t 



(6.10) 



with 



.4, 



l 2 



dr h(r, t)h n (r) 



+ 



dr h(r, t)h n (r) 



1/2 



(6.11) 



where the dot denotes the time derivative. To find the associated amplitudes, we consider 
the energy of a pure /i-excitation. Using (|6.5|) , and the boundary conditions on h, the 
quadratic energy is 

rL 

Integrating the second term by parts and using the equation of motion ( |6.7a| ) we find 



H h = f dr {{d t hf + {d r hf + AXh 2 } . 



(6.12) 



H h = 4tt J dr [{d t hf - h d 2 t h) 



(6.13a) 
(6.13b) 



Hence, the modulus squared of the amplitudes for this first mode are 
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|aiJ 2 = An A 2 uj ln 



(6.14) 



where u; ln = [4A + {x\ n / X) 2 ] 1 / 2 , x\ n = nn with n = 1, 2, • • • and the A n are given by ( |6.11|) . 
Equation ( |6.7b| ) produces an eigenmode whose non- vanishing components are of the form 

Vn(r,t) = y n (r) cosu 2n t, with 



Un{r) = A 2 „r ji(X 2n r) N 2n 



(6.15) 



where ui 2n = (1/2 + A^) 1 ^ 2 and A 2n = x 2n /L, with x 2n being the positive solutions to 
tan x 2n = x 2n (with this set of modes and those that follow, we will label the normal modes 
starting from n = 1). The parameters X 2n have been chosen so that y n (L,t) = 0, and the 
normalization constants N 2n are taken to be 

1/2 



2n 



Lsm x 2r 



(6.16) 



so that the y n (r) are orthonormal over [0, L\. To extract the amplitudes from a given solution 
we first expand the ^/-excitation as 



y(r,t) = J2 B nVn(r) cos u 2n t 



(6.17) 



with 



1 2 



dry(r,t)y n (r) 



+ 



^2n 



dry(r,t)y n (r) 



1/2 



(6.18) 



Using (|6.5| ), the quadratic energy of a pure y-excitation is 

H y = J o L dri [ (d t y) 2 + (d r y) 2 + ^ + 2 ^ . 
Integrating the second term by parts and using the equation of motion ( |6.7b|) we find 



(6.19) 



H v = 4n fdr {(d t y) 2 - y d 2 y} 



2 , .2 



UJ. 



2n 



(6.20a) 
(6.20b) 



Hence, the modulus squared of the amplitudes for the second mode are 
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\a>2n 



2 = knBluin 



(6.21] 



where cu 2n = [1/2 + [x 2n j X) 2 ] 1 / 2 , with x 2n being the positive solutions of tanx 2 „ = x 2n , and 
where the B n are given by ( |6.18| ). 

The remaining two modes are more involved since ( |6.7c| ) and ( |6.7d| ) are two coupled 
equations for if) and £. To disentangle these modes, we first rewrite these equations as 



d 2 ^-d^ + t + 2 4 + ^—,m + d t i] =0 



4 + r 2 



[d r i + dtl/j] = . 



r 2 r(4 + r 2 ) 

We now define ( = r(d r ib + dtC,)/ (4 + r 2 ), so that ( |6.22| ) may be rewritten as 



(6.22a) 
(6.22b) 



C , 2C 



dfC-d 2 C + i: + 



. 



(6.23a) 
(6.23b) 



Equation ( |6.23a| ) follows directly from ( |6.22a|) and the definition of (, while ( |6.23b| ) is derived 
as follows. First, take a time derivative of ( |6.22b|) . This gives a d 2 ip term in the square 
brackets, which may be eliminated using ( |6.23a|) to give 



2 ' r 2 r(4 + r 2 ) 



d r i + d r l[>' - 2C 



+ 4 = o, 



(6.24) 



where the dot and prime denote time and space derivatives respectively. We have written 
d 2 ^ rather than <9 t 3 £, d T ib' rather than d 2 ip, etc. for future convenience. Taking a spatial 
derivative of ( |6.23a| ) gives 



+ 



+ 2C' = . 



(6.25) 



Adding ( |B~24D and (|6~25| ), and using £ + — (4 + r 2 )(/r gives (|6.23b| ). 

These normal modes fall into two classes, one in which ( = and another in which 
( is non- vanishing. In the former case, ( |6.23a| ) may be solved for ib. We may then use 
dtC + d r ip = to solve for £. Thus, mode three takes the form ih 3n (r, t) = i/)3 n (r) sincj 3n t and 
£,3n(r,t) = (,3n(r) cosus n t, and after some algebra we find 
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ip3n(r) = A 3n r ji(X 3n r) N 3n 



A 



&n{r) = — [2ji(A 3n r) - A 3n r j 2 (A 3n r)] N 3n 

^3n 



(6.26a) 
(6.26b) 



where u 3n = (1/2 + Aj^) 1 ' 2 and A 3n = x 3n / L, with x 3n being the positive solutions to 
tanx 3n = x 3n /{l — x\ n ). The parameters A 3n have been chosen so that £ 3n (L,t) = (since 
£ vanishes, this automatically ensures that d r ip 3n (L,t) = 0). The normalization constants 
N 3n will be chosen below to ensure a convenient orthonormality relation for the ip 3n (r) and 
6n(r). 

We turn now to the other class of modes in which ( is non-vanishing. We can first solve 
( |6.23b ) for (, and then solve ( 6.23a ) for if) treating ( as a source. Then, using the definition of 
C, we can solve for £. Again, writing if>4 n (r, t) = if>^ n (r) smu^t and ^4 n (r, t) = ^4 n (r) cosu;4 n t, 
we find 



ip4n{r) = -TY- [2ji(A 4 „r) - A 4n r j (A 4 „r)] N 4n 

Ad 



A 4 n U) 'in 



-2X 4n r j (X 4n r) + 4(1 - \\ n )ji(\i n r) - 2\ An r j 2 {\i n r) 



N, 



in 



(6.27a) 
(6.27b) 



where a; 4n = (1/2 + A 2 n ) 1//2 and A 4ra = x 4n /L, with x 4n being the positive solutions 
to tanx4„ = X4„. The parameters A4„ have been chosen so that £,4 n {L,t) = and 
d r ip in (L,t) = 0, and the normalization constants N An will be chosen below. 
We expand a if>-^ excitation as 



J'=3,4 n 

£M) = C 'jnijnir) COS UJ jn t . 
i=3,4 n 



Choosing the normalization constants 

„2 \ „;„2 



A^4n 



2 x 3n) ^3n _|_ x 3 n X 3n 2 



2Lx 4 



3/i 



2Lx 2 



3n 



-1/2 



L 3 (2< + L 2 )sin 2 x 4 „ , L 5 (x 2 n -1) 



+ 



i" 



in 



x 



in 



-1/2 



the modes satisfy the orthonormality relations 



(6.28a) 
(6.28b) 



(6.29a) 
(6.29b) 
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rfr < - u jn £ jn (r) tOkm£,km{r) + — — 5- d r ip jn (r) d r tp km (r) 



10 



4 + r 



4 + r 



2 1 



L f 4 r 2 

dr <! uj jn ip jn (r) uJ km ip km (r) + — — - d r £ jn (r) d r £ k m(r) 



(6.30a) 



4 + r 2 



4 + r 2 



+ 2 On( r )Cfcm(^) 



(6.30b) 



Using ( p.30| ) in Q6.28|) , the overlap coefficients Cj n become 



Cj r 



drl — ijj 







jn 4 + r 2 



4 2 

+— — ^<9 r V>(r, i) d r ip jn (r) + — t/j(r, t)ip jn (r) 



4 + r 2 

L 



+ 



dr uj 







^ n 4 + r 2 



d t i>(r,t)il) jn (r) 



2\ 1/2 



(6.31) 



To extract the amplitudes, consider a pure excitation. Using (|6.5| ), the quadratic energy 
is given by 



e 2 + e 



+ 



4 + r 2 



^ 2 + 



2^ 2 £ 2 



J=3,4 n 



(6.32a) 
(6.32b) 



hence 



*jn I 



47T 



C 



J = 3,4 



(6.33) 



Even though we have solved the normal modes problem analytically in the continuum, 
the amplitudes \aj n \ 2 will be extracted using discrete numerical solutions. This is justified 
by the large size of our lattice: N = 2239, Ar = 0.04 (with A = 0.1). 

For computational purposes it is important to note that, strictly speaking, completeness 
sums involve all normal modes, but in a physically meaningful situation they will be satu- 
rated well before the normal mode indices reach the maximum value N. The highest normal 
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modes indeed correspond to artifacts of the discretization. Thus, to avoid unnecessary com- 
putational burdens, we will place a cutoff N mode ~ iV/5 to N/10 on the number of normal 
modes and calculate the Higgs and gauge boson particle numbers as 

■^V mode 

t'higgs = l°i«| 2 (6.34a) 

n=l 

^modc 

^gaugc = J2 \\ a 2n\ 2 + \a^n\ 2 + |a 4n | 2 } . (6.34b) 

n=l 

The total particle number is given by 

V = fhiggs + ^gauge • (6.35) 

We have verified that our results are insensitive to this cutoff, which means that short 
wavelength modes comparable to the lattice spacing are not excited in any appreciable 
manner. One should also note that our procedure for calculating the particle number is 
obviously gauge invariant (as it should be) since it makes use of an expansion into normal 
modes of gauge invariant variables. 
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FIG. 7. Decay of a small perturbation about the sphaleron: behavior of the particle number in the four 
normal modes of oscillation of the linearized system as function of time for lattice parameters ./V = 2239, 
Ar = 0.04 and N mo( ie = 200 with A = 0.1. The physical particle numbers are obtained by multiplying 
the asymptotic values in the graph by Air/g 2 ~ 30, which gives -/Vhig gs ~ 8 and -/Vg auge ~ 45, for a total 
physical particle number of iVph ys ~ 53. 
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In Fig. 7 we display the behavior of the particle number in the four normal modes 
of oscillation as function of time. The initial state is the small perturbation about the 
sphaleron in Fig. 5, which gives rise to outgoing spherical waves as the configuration decays. 
This is the state from which we start the stochastic sampling procedure described in the 
next section. Since the energy density is distributed over an expanding shell, the system 
quickly approaches the linear regime. This is apparent from Fig. 7 where, after an initial 
transition period in which the particle numbers of the four modes are not constant, they 
settle to values which are reasonably constant in time. We take this as evidence that the 
system has indeed reached an asymptotic linear regime where one can define a conserved 
particle number. 

There are two additional quantities that are useful in measuring the extent of linearity, 
namely the spectral energy e spec and the linearized energy ei in . The spectral energy is defined 
as the sum over normal mode energies, 

^modc 

^spec = X/ {^ln| a ln| + ^2n | | + ^3n|°3n.| +^4n|0-4n| } , (6.36) 
n=l 

while the linearized energy is defined by integrating the energy density in ( |6.5| ) expanded to 
second order in a perturbation about the vacuum, 

e lin = 4tt f°° dr{ (d t y) 2 + (d r y) 2 + % + t- + [d t hf + [d r hf + 4A/i 2 + 
Jo I r z 2 

jL- [ (ft£) . + mr ] + 1 + ^ [m-y- + (d^r] + ^} . 

Both the spectral and linear energy are gauge invariant since they have been defined us- 
ing gauge invariant quantities. If the system linearizes, then both e spec and eu n should be 
close to the conserved total energy e, which is given by the integral ( |6.5| ) (or in terms 
of gauge- variant variables by ( |3.7D ). The total energy of the configuration in Fig. 7 is 
given by e/An = 2.5447, while the asymptotic spectral and linear energies are given by 
Cspec/47r = 2.5679 and eu n /47t = 2.5685, and we see that the system has linearized to within 
one percent. (We also see that the sum over the energies of individual modes, although cut 
off at iV mo d e , essentially accounts for all the linearized energy.) 
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FIG. 8. Mode distribution for the asymptotic state of Figs. 5 and 7. This distribution is gauge invariant, 
and shows that all the particle are rather soft and comparable in energy. 

We can also investigate the mode distribution by examining the amplitudes \cij n \ 2 as a 
function of mode number n. As the system linearizes and the particle number becomes 
well defined, the mode distribution also becomes constant in time. Figure 8 illustrates the 
distribution of the asymptotic linear state of Figs. 5 and 7. Note that the population of the 
system is heavily weighted towards low lying modes. The mode cutoff used in calculating 
the particle number was iV mo de = 200, and we see that modes greater than about n = 150 
are not populated to any appreciable extent. The mode distributions are heavily peaked 
near n p k ~ 50, which corresponds to a frequency of co> p k ~ n v ^T\jL ~ 0.1. The perturbation 
about the sphaleron of Figs. 5 and 7 decays into about 50 rather soft particles (in physical 
units), each one of comparable energy. Finally, we point out that the mode distribution is 
gauge invariant as well. 

VII. STOCHASTIC SAMPLING OF INITIAL CONFIGURATIONS 

As we have discussed, our goal is to find the region in the e-v plane spanned by 
all topology changing classical solutions. More specifically, we would like to find the lower 
boundary of this region. The tools we have at our disposal allow us to vary the coefficients 
c n of ( |5.2|) , which defines the system as it passes over the sphaleron barrier, and to calculate 
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the corresponding energy e and incoming particle number v. From the computational point 
of view, e and v can be considered as known functions (albeit laboriously obtained) of the 
variables c n . We would then like to find 

Slower (e) = Min {cn;fixed e} v . (7.1) 

The particle number v may have several local minima since it is a highly non-linear 
function of the variables c n , and a straightforward constrained minimization procedure, such 
as a conjugate gradient technique, could fail to reveal the absolute minimum of v at a given 
e. We therefore decided to solve the problem using stochastic sampling. Stochastic sampling 
methods, driven by suitable weight functions and in combination with annealing techniques, 
have indeed proven very effective in exploring the overall structure of complicated surfaces 
and in approximating their global minima. 

Our procedure consists in generating "configurations" of the system weighted by a func- 
tion 

W = exp(-F) , (7.2) 

with 

F = (3e — jju . (7.3) 

By "configuration" simply we mean the collection of variables c n , which determine the whole 
evolution of the system. Since e and v are functions of c„, the weight given by ( |7.2| ) and 



( [7.31 ) is also a function of c n and defines a probability distribution 

dP = Z- l l[dc n W(c n ) . (7.4) 

n 

We will generate topology changing configurations distributed according to (|7.4|). Clearly, by 
taking large values for the parameters f3 and fi we will drive the distribution strongly towards 
the lower boundary in the space of all topology changing solutions. By using different ratios 
fi/(3 we will be able to drive the distribution in different direction and thus follow the lower 
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envelope of the region, while temporarily lowering the values of /i and/or (3 will allow us to 
anneal the distribution. We will typically take f3 between 50 and 1, 000 while /i will range 
between 1, 000 and 20, 000. 

To generate the desired distribution we have used a Metropolis Monte Carlo algorithm. 
Starting from a definite configuration c n , we randomly select one of the variables c, and 
perform a variation q — > = q + Aq (in our computation, the Aq are Gaussian distributed 
with a mean of 0.0008.) The system is evolved backward and forward in time and we 
calculate the energy, in-state particle number and change of winding number. If the winding 
number does not change, we proceed to vary another of the variables c n . If the topology 
changes, we evaluate AF = j3Ae + fiAu and the new value d i is accepted with conditional 
probability p = Min[l, exp(— AF)]. Specifically, we generate a pseudorandom number r 
uniformly distributed between and 1, and if r < exp(— AF) the change is accepted and 
the new value replaces the old one. Otherwise, if r > exp(— AF) the old value is kept and 
we select another of the variables c n for a possible upgrade. (We should note here that when 
the winding number changes, even if the trial value d { is rejected, we still record its value 
and the corresponding values of e and v, since they do correspond to a possible topology 
changing evolution.) 

It must be emphasized that although our algorithm generates a distribution of topology 
changing solutions of the equations of motion, this distribution represents only a computa- 
tional device and carries no special physical significance. Indeed, the probability measure 
( [7.4| ) is based on the arbitrary choice of variables c n and no Jacobian factor of any kind 
has been introduced. It would be possible to define a measure which represents a phys- 
ically meaningful distribution, and our notation f3 and \x for the weights of e and v has 
been inspired by the analogy with a grand canonical ensemble. But still, in the present 
context, there is no reason for defining any particular physically meaningful measure and no 
justification for the attached computational costs. 
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FIG. 9. Monte Carlo results with lattice parameters of TV = 2239, Ar = 0.04 (giving L=89.56), 
TVmode = 200 and iV sp h = 50, and with a Higgs self-coupling of A = 0.1. The solid line marks the 
sphaleron energy e sp h = 47r(2.5426), below which no topology changing process can lie. The diamond 
represents the configuration from which we seeded our Monte Carlo search. To obtain quantities in physical 
units, multiply the numbers along the axes by Air/g 2 ~ 30. The energy axis extends from about 10 TeV 
to 15 TeV, while the particle number axis ranges from about 30 particles to 60. 

Figure 9 illustrates the results of our Monte Carlo investigation. It represents about 300 
hours of CPU time on a 16 node partition of a CM-5. We generated approximately 30,000 
configurations of which approximately 3,000 representatives are plotted in the figure. We 
have chosen lattice parameters TV = 2239 and Ar = 0.04, for a lattice of size L = 89.56. 
We have used a cutoff iV mo de = 200 on the sums over the modes, and the dimension of the 
initial configuration space over which we have sampled is determined by yV sph = 50. We 
have taken the Higgs self-coupling to be A = 0.1, which in lattice units corresponds to a 
mass of about m# = (40Ar) _1 , or a physical mass of run = 110 GeV. As one can see, our 
lattice is sufficiently dense that there are many lattice sites within a single Higgs Compton 
wavelength. 

It is apparent from Fig. 9 that our search procedure is effective in reducing the particle 
number and in exploring the lower boundary of the space of topology changing classical 
solutions. The complex nature of this space (or at least of our search procedure) is also 
apparent from the figure, in that one can clearly observe two breaks in the outline of the 
lower boundary at e/47r ~ 3 and e/4n ~ 3.4. The reason for the discontinuity is that in 
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a first extended search we did not verify that every individual solution changed topology 
(performing this check is costly in computer time), trusting that topology change would be 
the typical outcome of an evolution which passes over the sphaleron barrier. A subsequent 
analysis revealed however that for a whole subset of our configurations, comprised between 
e/47r ~ 3 and e/47r ~ 3.4, the topology did not change: the system went over the sphaleron 
barrier a second time in the reversed direction and returned to the original topological sector. 
We discarded all these configurations and verified that the topology changed in all the 
remaining ones. We then implemented the check for topology change at every Monte Carlo 
step and restarted our sampling procedure by annealing a topology changing configuration 
obtained for e/47r ~ 3. This second search produced the set of configurations which stand 
out at slightly lower v between e/47r ~ 3 and e/47r ~ 3.4. 



FIG. 10. Topology changing transition obtained after many Monte Carlo iterations: behavior of the x 



Our search procedure not only leads to classical solutions with lower particle number, 
but is effective in selecting configurations with special properties in the in-state (these two 
features of course go hand in hand). In Fig. 10 we illustrate the entire evolution for one 
of the topology changing processes with low particle number, corresponding to one of the 
points at the bottom-right corner of the plot in Fig. 9. Figure 10 should be contrasted 




field. For a full color figure see http: / /cthulu.bu.cdu/~bobs /bviolatc.htm] . 
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with Fig. 6 in which the evolution of our Monte Carlo seed configuration is illustrated. 
The change is dramatic. The in- and out-states in Fig. 6 look rather symmetric, up to the 
topology change of the out-state concentrated about the origin. The particle numbers of the 
in- and out-states are about the same and of order 50 in physical units [yijAix = 1.747 and 
v /Air = 1.750, respectively). Figure 9 shows that after many Monte Carlo iterations we have 
managed to filter initial configurations c n so that the in-state particle numbers are about 
40% lower (ui/An ~ 1.10), and from Fig. 10 it is apparent that the in-states are now much 
different from the out-states. The former are narrow with the spectrum shifted towards 
shorter wavelengths, while the outgoing states still display the broad long-range waves seen 
at both ends of the evolution in Fig. 6. Indeed, the particle number in the out-state remains 
high. 
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FIG. 11. Behavior of the in-state particle number in the four normal modes. The initial state was 
obtained after many Monte Carlo iterations, and soon linearizes. Note, however, that mode 3 remains about 
10% non-linear. 

More details of the configurations selected by our sampling procedure are revealed by 
Figs. 11 and 12. Figure 11 illustrates the behavior of the particle numbers associated with the 
four normal modes for the initial configuration c n used to generate Fig. 10. The asymptotic 
particle numbers in Fig. 11 are associated with the in-state of the physically relevant time 
reversed solution of Fig. 10. Figure 12 illustrates the mode distribution of this in-state. 
These figures should be contrasted with Figs. 7 and 8 which display the same quantities 
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at the beginning of our search. The change is again very impressive. In particular, it is 
clear that the stochastic sampling procedure has selected classical solutions where the mode 
distribution in the in-states is shifted towards higher frequencies and shorter wavelengths. 
Of course, this is necessary for a reduction of the ratio vje. 
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FIG. 12. Mode distribution of the asymptotic in-state of Fig. 11. 

Although our results show a marked decrease in the particle number of the incoming 
state, nowhere in the energy range we have explored does v drop below 47r, or in physical 
units iVp hys > 30 for E < 15 TeV. This is a far cry from the value Ap h y S = 2 which would be 
needed to argue that baryon number violation can occur in a high energy collisions. From this 
point of view our present results are limited and should be pushed to much higher values of 
e. In the next section we will make some comments about our future plans to explore higher 
energies and discuss other investigations which can shed further light on the properties of 
the system. As of now the computational resources at our disposal, together with the rather 
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ambitious number of points we have used for our numerical study, have not permitted us to 
go beyond the energy range we have explored. We believe that our results, as well as the 
formalism we have established, are nevertheless interesting enough to warrant publication. 
In some respect, the choice of a number of points as large as our current N = 2239 has 
been an error of strategy. In a preliminary investigation, described in [17], we had used 



N = 256. The number of points in the lattice determines of course the ultraviolet cutoff and 
this in turn implies a minimum value for the ratio vjt. This quantity is indeed minimized 
by placing all the weight in the highest mode iV mo d e , giving (V/e) m m = 1/^modc ~ L/N^o^. 
With iV = 256 points we saw the onset of this constraint, and we decided to chose a lattice 
size that would push the lower limit on vjt to a much smaller value closer to the physically 
relevant domain. With the parameters of our present calculation, the minimum would occur 
at (^/e) m i n ~ 0.15. However, the increased computational burden, together with the fact 
that the stochastic sampling moved in the e-v plane at a much slower rate than we had 
anticipated, prevented us from saturating this lower bound. 

It is still interesting to extrapolate our results to obtain information about the possible 
behavior of the boundary in the e-v plane of topology changing solutions. For this purpose 
we binned all our data into subintervals of width Ae = 0.005. Within every bin we selected 
the point with lowest v. We then fitted these points to the hyperbola 

{y - ae - c x ){v - v^) = c 2 , (7.5) 

where a and are the free parameters of the fit. The quantities C\ and c 2 (which are 
constants with respect to v and e but depend on a and v^) are given by C\ = 2z/ sph — ae sp h— Voo 
and C2 = — (z/ sp h — i^o) 2 , where e sp h = 2.5447 and z/ sp h = 1.7478 are the energy and particle 
number in the limiting case in which the configuration approaches the sphaleron itself (in 
practice the e and v of the configuration from which we started the Monte Carlo search). 

This fit is motivated by simple physical considerations. We would expect the lower 
boundary of the region of topology changing transitions to saturate either at v = or at 
some finite value of v. The boundary of the domain must go through the sphaleron and 
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should have an infinite slope there. Indeed, since the topology changing classical solutions 
become complex when e decreases below e sp h, one would expect the boundary curve v = v(e) 
to have a square root singularity at e = e sp h- Finally, although the upper boundary of the 
region is of little interest to us, it is not unreasonable to parameterize it in terms of a straight 
line of constant slope. This is the line one would find if the upper bound were obtained by 
putting all the energy in a single mode of frequency to (in which case the slope a = 1/uj), 
or since this is unrealistic, if the the mode distribution could be well approximated in terms 
of some effective frequency u e s = 1/a. The hyperbola of (|7.5| ) is the simplest curve with all 
these properties. 

The results of our fits are shown in Figs. 13 and 14. In Fig. 13 all the data points 
have been used, and the solid line represents the unconstrained fit while the dashed line is 
obtained by requiring = 0. Since one can argue that what ought to be fit is the lower 
boundary of the region, and that insofar as our points display a slight discontinuity and 
cannot all belong to this boundary, we have repeated the fit removing all the points which 
lie above the unconstrained fit in Fig. 13. The results of this second fit are reproduced in 
Fig. 14. For Fig. 13, the unconstrained fit has parameters a = 0.257 and = —0.294 while 
the constrained fit has a = 0.319. Figure 14 has the parameters a = 0.238, = —0.530 
and a = 0.341 respectively for the unconstrained and constrained fits. 




2 4 6 8 10 

e (units of 4ji) 

FIG. 13. Hyperbolic fits to full data set. The asymptotic particle number is constrained to vanish for 
the dashed line, while it remains unconstrained for the solid line. 
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FIG. 14. Same as Fig. 13, except the data set was reduced by those point lying above the previous 
unconstrained fit. 

It is interesting to observe that that the unconstrained fits lead to an asymptotic value for v 
smaller than zero, which shows that one cannot read any indication of a lower bound on the 
particle number in our present data. Our results cover a range of energies which is too small 
to derive any reliable conclusion about whether and when the particle number could reach 
the value two. One can nevertheless insert physical units in the results of our fits and see at 
what energy values the incident particle number would become equal to two. This simple 
exercise gives energies of 110.37 TeV and 447.20 TeV respectively for the unconstrained and 
constrained fits of Fig. 13, and energies of 75.06 TeV and 418.61 TeV for the corresponding 
fits of Fig. 14. 

We conclude this section with a few technical remarks. Since our entire procedure is 
based on the calculation of the particle number after the system has reached the linear 
regime, we should make sure this quantity is evaluated in a reliable manner. Now, it is 
clear from the graphs of Figs. 7 and 11 that, while the particle number becomes reasonably 
constant towards the end of the evolution, it still exhibits oscillations possibly as large as 
10%. This might cast doubts on the validity of our stochastic sampling technique, where the 
steps in initial parameter space induce variations of the particle number as small as 10~ 4 . 
The solution we have adopted is to define a "computational particle number" v c (which is the 
quantity represented in Fig. 9). With a lattice of infinite spatial extent, even in presence of an 

53 



ultraviolet cutoff arising from finite lattice spacing, and barring the existence of conservation 
laws giving rise to particle phenomena, the system will eventually linearize fully and the true 
particle number v will be well defined and constant to any degree of precision. Since we begin 
with an initial state localized around the origin, we may conceptually think of this as being 
defined over an infinite lattice, although in practice we use a lattice of finite extent. Thus, 
every initial configuration c n conceptually determines a unique particle number v. This may 
not be accessible to us, but it exists. We define a quantity v c which we can measure as 
follows: we evolve the system for a definite amount of time T and then for an additional 
time AT (in our calculation To = 60 and AT = 8). Over the interval To, Td + AT we measure 
the particle number at times T = T . . . T m chosen at random (in our calculation we take 
m = 10 and choose T x through T 10 to be 61.55, 62.51, 63.27, 63.70, 64.77, 65.25, 65.33, 
65.71, 66.59, 68.00 respectively), but fixed for the entire calculation. The computational 
particle number v c is defined as the average of the particle numbers measured at Tj. Again, 
v c is a well defined function of the parameters c n , and uniquely determined by this initial 
configuration. The crucial point is that v c tracks v. The quantities v and v c may differ 
by as much as 10%; however, if we reduce v c by a certain factor, we can be confident that 
the true particle number v has also been reduced by the same factor, up to a relative error 
given by the approximation by which v c tracks v. Finally, we should make sure that v c 
is, computationally, a well behaved function of the parameters c n , i.e. that the functional 
relation between the chosen c n and the measured value of v c is not spoiled by numerical 
errors. This we have verified explicitly. On a sample configuration we have stepped every 
individual parameter c n by values an order of magnitude smaller than the typical steps in 
our stochastic sampling procedure and have verified that the corresponding changes in v c 
are regular and well accounted for by the first few terms of a Taylor series expansion in Ac n . 

VIII. CONCLUSIONS 

We have developed a computational procedure that allows us to explore the space 
of classically allowed topology changing transitions leading to baryon number violation. 
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With our method we have been able to trace the lower boundary of the region spanned by 
topology changing evolution in the energy versus incoming particle number plane, up to 
energies approximately one and a half times the sphaleron energy and with a reduction of 
the incoming particle number by approximately 40%. The corresponding solutions display 
dramatically different features in their incoming state from the solution used to seed the 
Monte Carlo search (in which there was just barely enough energy to cross the sphaleron 
barrier), one of the most notable differences being a marked shift of the in-state spectral 
mode distribution towards higher frequencies and shorter wavelengths. Within the domain 
we have explored there is no indication of an emergent lower limit on the particle number of 
the incoming state. Indeed, a hyperbola fit to our data, motivated by the expected physical 
properties of the boundary of the domain of topology changing evolution, is quite compatible 
with a zero lower bound on the incoming particle number. 

Our results are unfortunately rather limited in the extent of energy and particle num- 
ber which we have been able to explore. However impressive may be the change in the 
properties of the solutions spanned by our search, the fact remains that the lowest particle 
number we have been able to reach is, in physical units, approximately 30. An even more 
serious shortcoming of our results is that our method can only establish an upper bound on 
the minimum particle number at any given energy: when our search produces a topology 
changing solution of given e and v>, it establishes by construction that the lower boundary of 
the classically allowed transitions cannot lie above that point, but we cannot rule out that 
it might lie substantially below and that the stochastic search simply failed to come close 
to it. 

However, the mere fact that the analytically intractable non-linear equations of motion 
are amenable to a reliable computational solution is, we believe, a very important result, 
perhaps the most important fact emerging from our analysis. By solution, we mean much 
more than just the implementation of a numerical integration algorithm of the evolution 
equations. Our study makes it clear that a whole range of detailed questions about the 
entire space of solutions can be tackled and solved by computational means. 
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The results we have established thus far naturally lead to further investigation. By in- 
vesting more computational resources it will be straightforward to extend the exploration to 
substantially larger energies. However, one can do more than that. The detailed information 
obtained about the spectral composition of the incoming states with low particle number 
suggests that that one may explore the properties of such states directly. For instance, one 
could try to shift the mode distribution further towards shorter range, while verifying that 
the ensuing evolution still changes topology. This runs somehow against our original notion 
that it would be very difficult to start from the selection of the incoming state and still ob- 
tain a topology changing solution, but now we are no longer dealing with a blind sampling of 
incident states. From this point of view we find very inspiring some recent results obtained 
by Farhi, Goldstone, Lue and Rajagopal who, in a study of collision induced soliton decay, 
were able to produce the "unwinding" of the soliton and its subsequent decay by directing 
against it waves which carry a short range twist of the phase of the complex field (we refer 
to the original work of Ref. |L8] for an elucidation of this possibly cryptic sentence). It is 
interesting that in computer animation which we generated to clarify the properties of the 
evolution, we have seen analogous twists in the phase of the 2-dimensional field x m the 
asymptotic states. Of course one must be careful in defining effects which pertain to gauge 
variant quantities, but a careful study of the properties of the asymptotic states may provide 
important clues for understanding the mechanisms leading to classically allowed transitions 
with low incoming particle number. 

Finally, a complementary approach to the study of classically allowed transitions consists 
in studying the classically forbidden processes. As we have already mentioned in the Intro- 
duction, a very powerful formalism for the study of such processes has been established in 



Ref. and applied recently in Ref. |19| to the study of collision induced decay of the false 
vacuum. The method of Ref. || requires that one solves analytically continued equations 
of motion along a suitable contour in the complex time plane and that one imposes bound- 
ary conditions based on the normal mode expansion of the fields in the linearized domain. 
Thus a large part of the formalism we have developed in this paper will carry over to the 

56 



study of classically forbidden processes, and we plan to make this the subject of a future 
investigation. 
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